Bayesian Statistics

Loading required package: plotly
Loading required package: ggplot2
Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

Loading required package: scales

A Bayesian is one who, vaguely expecting a horse, and catching a glimpse of a donkey, strongly believes he has seen a mule

Why might someone use Bayesian rather than traditional frequentist?

• If they have prior beliefs as to the truth before an experiment is conducted
• If you want to predict how the distribution of a population might change, rather than a point estimate. For example, traditionally, we might try and predict the mean number of people who choose to order after going on our website.

Quick recap of the ‘Classical Approach’ (‘Frequentist’)

set.seed(0)
p <- 0.5
n <- 10

Let’s say we have a coin, and we want to know the probability that we flip heads. Frequentists assume that there is only one true, fixed, unknown probability that the coin yields heads. For ease of notation, let’s call the constant parameter \(p(\text{head}) = \theta\).

To try and discover the true probability of getting a head, we flip the coin 10 times, and count the number of heads to estimate the probability. Frequentists believe that if we repeated this sampling proceedure of 10 flips many, many times, we should expect the average across all of our samples to be reflective of the true population value (hence why it is called ‘frequentist’).

# Create repeated samples
num_repeated_samples <- 10^3
repeated_samples <- rbinom(n = num_repeated_samples, size = n, prob = p)
repeated_samples.cum_p_mean <- cumsum(repeated_samples)/cumsum(rep(n, num_repeated_samples))

# Cumulative mean tends to p
plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = 1:length(repeated_samples), y = repeated_samples.cum_p_mean, name = 'Cumulative mean with sampling') %>%
  layout(showlegend = F
         , yaxis = list(range = c(0, 1))
         , title = "Cumulative sample means tends to true p") # , 

k <- round(n*p*0.5,0) # make up fake k that is 'unlikely'
sample_mean <- k/n

binomial_probability <- function(k, n, p) {
  n_choose_k <- factorial(n)/(factorial(k)*factorial(n-k))
  # NOTE - equivalent to gamma_function(k+n-k)/(gamma_function(k)*gamma_function(n-k), or 1/B(k,n-k)
  prob_sample <- n_choose_k*(p)^k*(1-p)^(n-k)
  return(prob_sample)
}

k_heads_maxprob <- round(binomial_probability(k,n,k/n),4)*100
k_heads_actprob <- round(binomial_probability(k,n,p),4)*100
np_heads_actprob <- round(binomial_probability(n*p,n,p),4)*100

But in the real world, it’s not possible to take lots of samples over and over again - and we generally just take one. Taking a sample one time, we might observe 2 heads. This is equivalent to creating a sample of size \(10\) and observing a mean probability of \(0.2\).

We don’t know what the true probability of flipping heads is (\(\theta\)), but we can plot the likelihood of observing 2 heads for every possible \(\theta\) the coin might be. On the left graph below, for every possible value that the probability of flipping heads could be on the x axis (somewhere between 0 and 1), the likelihood of observing \(2\) heads is plotted on the y axis, giving the orange line. For example, the likelihood of observing \(2\) heads if the probablity of flipping heads is \(0.2\) is \(30.2\%`\). If the probablity of flipping heads is actually \(0.125\), then the likelihood of flipping 2 heads in our sample is only \(4.39\%\).

It appears most probable from our sample then that our true probability is indeed 0.2 - the most likely probability is what we observe as the mean in our sample (the maximum likelihood estimation is rsample_mean). Thus, a frequentist would conclude that the most likely value of \(\theta\) is \(0.2\), and that the true value of \(\theta\) is hopefully near to it.

However, despite treating the true \(\theta\) as a fixed (non-random) value, frequentists acknowledge there is random sampling error. Even if the probability of flipping heads remaining constant, sometimes there will be 2 heads, sometimes there will be 10 heads etc simply due to chance. We can illustrate this by plotting the likelihood of observing different numbers of heads given a fixed value of \(\theta\). In the graph on the right, the probability of getting different numbers of heads between 0 and 10 is plotted in green if the value of theta is \(0.2\), and blue if the probability is \(0.125\).

In this way then, despite believing the actual (population) probability \(\theta\) to be fixed at one value (and from the sample, the most probable being \(0.2\)), the graph on the right illustrates how frequentists believe that the number of heads observed is a random variable (conditioned by \(\theta\)). By random variable, this means that the number of heads could be many different values before sampling (0-10) rather than being pre-determined/fixed as how they treat \(\theta\). The likelihood of observing any of these possible values is conditioned by the fixed probability \(\theta\) as per the distribution on the right (if the true value of \(\theta\) is \(0.125\), then observing \(1.25\) heads in our sample is a lot more likely than observing \(10\), and observing \(2\) heads is very unlikely).

# for all the different values of theta, p(H=2|θ)
theta_prob_k_heads <- sapply(seq(0,1,0.01), FUN = function(p) binomial_probability(k=k,n=n,p=p))
# for all the different values of outcome k, p(H=k|θ=0.2)
likelihood_est <- sapply(seq(0,10,1), FUN = function(s) {binomial_probability(k=s,n=n,p=k/n)})
# for all the different values of outcome k, p(H=k|θ=0.5)
likelihood_act <- sapply(seq(0,10,1), FUN = function(s) {binomial_probability(k=s,n=n,p=p)})


subplot(
  plot_ly(type = 'scatter', mode = 'lines') %>%
    add_trace(x = seq(0,1,0.01), y = theta_prob_k_heads, name = paste0('P(H=',k,'|θ)')) %>%
    layout(xaxis = list(title = 'θ'))
  , plot_ly(type = 'scatter', mode = 'lines+markers') %>%
    add_trace(x = seq(0,10,1), y = likelihood_est, name = paste0('p(H=k|θ=',k/n,')')) %>%
    add_trace(x = seq(0,10,1), y = likelihood_act, name = paste0('p(H=k|θ=',p,')')) %>%
    layout(xaxis = list(title = '# heads'))
  , shareY = T
  , nrows = 1)

# 
# plot_ly(type = 'scatter', mode = 'lines') %>%
#   add_trace(x = x, y = likelihood_est, name = 'Estimated from sample') %>%
#   add_trace(x = x, y = likelihood_act, name = 'Actual likelihood') %>%
#   layout(title = "Likelihood function: X ~ Binom(n,p)"
#          , xaxis = list(title = 'theta')
#          , yaxis = list(title = 'likelihood'))

Hence, to acknowledge that there is randomness associated with this sample mean, frequentists construct confidence intervals. Again to follow the frequentist interpretation, a 95% confidence interval means that, if the sampling proceedure was repeated many, many times, that 95% of the time the true value of \(\theta\) would fall within the confidence interval. (Nuanced note - this is not the same as saying that, in one sample, there is a 95% chance that the true value of \(\theta\) lies within it. This is because \(\theta\) is a fixed value, not a random variable, so it either falls within the sample confidence interval or not.)

true.p.se <- sqrt(p*(1-p)/n)

## BINOMIAL (NORMAL APPROXIMATION)
# The larger the number of successes, np, and failures, n(1-p), the better the normal approximation for confidence intervals (p near 0.5, n as large as possible)
sample.binom.k <- rbinom(n = num_repeated_samples, size = n, prob = p)
sample.binom.p <- sample.binom.k/n
  ## population standard error/confidence interval
    true.binom.p.conf.02.5 <- sample.binom.p-1.96*true.p.se
    true.binom.p.conf.97.5 <- sample.binom.p+1.96*true.p.se
    true.binom.theta_in_conf <- true.binom.p.conf.02.5 < p & true.binom.p.conf.97.5 > p
  ## sample standard error/confidence interval
    sample.binom.p.se <- sqrt(sample.binom.p*(1-sample.binom.p)/n)
    sample.binom.p.conf.02.5 <- sample.binom.p-1.96*sample.binom.p.se
    sample.binom.p.conf.97.5 <- sample.binom.p+1.96*sample.binom.p.se
    sample.binom.theta_in_conf <- sample.binom.p.conf.02.5 < p & sample.binom.p.conf.97.5 > p

## BETA DISTRIBUTION (CONTINUOUS VERSION OF BINOMIAL)
# Since the binomial distribution is non-continuous, we can get confidence intervals closer to 95% from the equivalent beta distribution for the same np & n(1-p)
sample.beta.k <- rbeta(n = num_repeated_samples, shape1 = p*n, shape2 = (1-p)*n)*n
sample.beta.p <- sample.beta.k/n
  ## population standard error/confidence interval
    true.beta.p.conf.02.5 <- sample.beta.p-1.96*true.p.se
    true.beta.p.conf.97.5 <- sample.beta.p+1.96*true.p.se
    true.beta.theta_in_conf <- true.beta.p.conf.02.5 < p & true.beta.p.conf.97.5 > p
  ## sample standard error/confidence interval
    sample.beta.p.se <- sqrt(sample.beta.p*(1-sample.beta.p)/n)
    sample.beta.p.conf.02.5 <- sample.beta.p-1.96*sample.beta.p.se
    sample.beta.p.conf.97.5 <- sample.beta.p+1.96*sample.beta.p.se
    sample.beta.theta_in_conf <- sample.beta.p.conf.02.5 < p & sample.beta.p.conf.97.5 > p

# First one-time sample confidence intervals
one_time.sample.p.se <- sqrt((k/n)*(1-k/n)/n)
one_time_sample.conf_02.5 <- (k/n)-1.96*one_time.sample.p.se
one_time_sample.conf_97.5 <- (k/n)+1.96*one_time.sample.p.se

# One time sample either contains true value of theta or not
plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = rep(c(one_time_sample.conf_02.5, one_time_sample.conf_97.5), each = 3), y = c(1,-1,0,0,-1,1), name = 'Confidence interval') %>%
  add_trace(x = p, y = 0, mode = 'markers', name = 'True value', marker = list(size=50, symbol = 'x')) %>%
  layout(showlegend = F
         , yaxis = list(range = c(-2, 2), zeroline = F, showticklabels = F)
         , xaxis = list(range = c(one_time_sample.conf_02.5 - 0.1, one_time_sample.conf_97.5 + 0.1))
         , title = "One-time sample confidence interval either contains true value of theta or not")

# Cumulative confidence intervals contain 95% of the time
plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = 1:num_repeated_samples
            , y = cumsum(true.beta.theta_in_conf)/cumsum(rep(1,num_repeated_samples))
            , name = 'Cumulative mean with sampling') %>%
  layout(showlegend = F
         #, yaxis = list(range = c(0, 1))
         , title = "Cumulative confidence intervals tends to contain p 95% of the time")

NA

In order to obtain these results, certain conditions have to be met, including:

• The randomly sampled observations are independent, so observing any one value does not have an impact on what subsequent values we observe might be. In our example, flipping a head once does not impact the likelihood of flipping it again. (If running an experiment in practice, this might mean that our population is sufficiently large when compared to our sample that we don’t need to be concerned with replacement. Another thing to watch out for is network effects - imagine an airline increases price of seats as they sell. Say we want to experiment in offering some customers a special offer for airline tickets, our variant, and some not, our control: the customers with the special offer might be more likely to buy, and this will shorten the supply, increasing the price for those without the special offer, and decreasing their likelihood to buy).
• Any sampled observations are identically distributed to the population (e.g. if we take our sample from the USA, they might have a different propensity to watch TV than the UK.)
• Because it is a random sample, the values observed are characterised by a probability distribution associated with \(\theta\). (e.g. we don’t know what the number of heads oberved will be before sampling, but the likelihood of what we observe is pre-determined through its associated probability distribution with \(\theta\))

Okay, so what exactly does a Bayesian approach mean?

Bayesian inference is broadly similar to frequentist inference. The key difference though is that, under a Bayesian approach, the true value of \(\theta\) is treated as a random variable, rather than a fixed constant.

So, for example, we might believe that in our economy, \(90\)% of coins are fair, and \(10\)% of coins are biased, with biased coins flipping heads \(20\)% of the time.

Under the guise of being frequentists previously, we had no prior belief as to what \(\theta\) would be: we only made an educated guess based on the 2 heads we observed. We believed that there is one fixed value of the true \(\theta\), and that there was only randomness associated with sampling, so we created a confidence interval to portray this.

Now, as bayesians, we are also incorporating our prior beliefs as to what \(\theta\) is - that it is either fair or biased. We believe that there is not just random error associated in sampling - we also treat \(\theta\) as an additional random variable, since we don’t know whether our coin is fair or biased either, and fair coins will have a different sample mean distribution than biased ones.

fair_coin_k_prob <- binomial_probability(k, n, p)
biased_coin_k_prob <- binomial_probability(k, n, k/n)

Given that we observe 2 heads, the likelihood of this occuring if the coin is fair is \(4.39%\)%, and \(30.2%\)% if the coin is biased. In the absence of informative prior information then, it might appear that the coin is more likely to be biased.

However, we also need to take into account that \(90\)% of coins in the population are fair, so even before observing 2 heads, we have some belief about whether the coin is fair or biased (i.e. the probability distribution associated with the random variable \(\theta\)).

fair_coin_prob <- (fair_coin_k_prob * 0.9)/((fair_coin_k_prob * 0.9) + (biased_coin_k_prob * 0.1))
biased_coin_prob <- (biased_coin_k_prob * 0.1)/((fair_coin_k_prob * 0.9) + (biased_coin_k_prob * 0.1))

# data.frame(
# matrix(data = c(binomial_probability(2,10,0.5)*0.9 # p(fair, 2 heads)
#                 , (1-binomial_probability(2,10,0.5))*0.9 # p(fair, not 2 heads)
#                 , 0.9 # p(fair)
#                 , binomial_probability(2,10,0.2)*0.1 # biased, 2 heads
#                 , (1-binomial_probability(2,10,0.2))*0.1 # biased, not 2 heads
#                 , 0.1
#                 , binomial_probability(2,10,0.5)*0.9 + binomial_probability(2,10,0.2)*0.1
#                 , (1-binomial_probability(2,10,0.5))*0.9 + (1-binomial_probability(2,10,0.2))*0.1
#                 , 1
#                 )
#        , nrow = 3
#        , dimnames = list(c('Flip 2 heads', 'Not flip 2 heads', 'Total')
#                          , c('Fair', 'Biased', 'Total'))
#        )
# )

..and so we can find the likelihood we got a fair coin, given the fact we observed 2 heads:

\[ P(\text{fair}|H=2) =\frac{P(\text{H=2 and fair})}{P(\text{H=2 and fair})+P(\text{H=2 and biased})} =\frac{P(\text{H=2|fair}) \times P(\text{fair})}{P(\text{H=2})} \]

\[ P(\text{fair}|H=2) =\frac{0.0439 \times 0.9}{0.0439 + 0.302}=0.567 \]

\[ P(\text{biased}|H=2) =\frac{P(\text{H=2 and biased})}{P(\text{H=2 and fair})+P(\text{H=2 and biased})} =\frac{P(\text{H=2|biased}) \times P(\text{biased})}{P(\text{H=2})} \] \[ P(\text{biased}|H=2) =\frac{0.302 \times 0.1}{0.0439 + 0.302}=0.433 \]

In other words, our prior belief (before sampling and observing any data) was that the likelihood of our coin having \(\theta\) equal to \(0.125\) was \(90\%\), and the likelihood of it having theta equal to \(0.2\) was \(10\%\). After observing 2 heads, we update our beliefs, so that the likelihood of the probability of our coin flipping heads (\(\theta\)) being equal to \(0.5\) is now equal to \(56.7\%\), and of it being equal to \(0.2\) being \(43.3\%\)

It is this incorporation of prior data that is either seen as Bayesian’s biggest advantage or pitfall. On the one hand, experiments are not abstract devices, and some knowledge about the process being investigated before obtaining the data is known and arguably should be incorporated. On the other, incorporating subjective opinions, particularly strong ones, may mean that you do not learn the true values you are trying to derive. Bayesian is thus analysis that uses a set of observations to change opinion rather than as a means to determine ultimate truth.

To help describe bayesian analysis, the following terms are often used:

• Prior distribution: \(Pr(\theta)\) - represents existing belief about \(\theta\) (Represents what was thought before seeing the data). For example, in the binomial coin example above, we have prior knowledge that 90% of coins in our economy are fair, so \(P(\theta)=P(Biased)=0.1\)
• Evidence: \(X\) - what we just observed (\(2\) heads)
• Marginal probability: \(Pr(X)\) - the total probability of the data across all possible values of the parameter \(\theta\). It doesn’t actually depend on \(\theta\) and isoften referred to as the proportionality factor/normalising constant (it makes sure all the scenarios are modelled. For example, \(P(H=2) = P(H=2|Fair) \times P(Fair) + P(H=2|Biased) \times P(Biased)\))
• Likelihood function: \(Pr(X|\theta)\) - the probability of \(\theta\) given the data observed (Represents the new data available). For example, \(P(H=2|Biased)=0.3019899\)
• Joint probability density function: \(Pr(X,\theta)=Pr((X|\theta).Pr(\theta)\) - used to modify prior beliefs through Bayes Theorem
• Posterior distribution: \(Pr(\theta|X)\) - the posterior density represents the knowledge about the model parameters after observing the data (Represents what is now thought given both prior data and data just obtained)
• Conjugancy: occurs when the posterior distribution is in the same family of probability density functions as the prior belief, but with new parameter values (updated to reflect what has been learned from the data). The posterior comes from the same family as the prior, rather than the data’s distribution.

And thus, we can then describe Bayes Theorem:

\[ P(\theta|X) =\frac{P(X|\theta) \times P(\theta)}{P(X|\theta) \times P(\theta) + P(X|\theta') \times P(\theta')} =\frac{P(X|\theta) \times P(\theta)}{P(X)} \]

Okay - the example above might not quite be a fair comparison, since the Bayesian example had quite a strong prior: that there were only two possible values for theta, and that 90% of the time it is 0.5.

A more ‘objective’ prior that is a fairer comparison to the frequentist example is to give all values of theta have an equal likelihood - in other words, a uniform distribution. We can model this using a beta distribution. Going forwards too, let’s make this more generic, and term the number of heads observed in the sample as \(k\).

A quick introduction to the beta distribution: it can model lots of different functional forms using \(a\) and \(b\) as parameters, by inputting them into the following function:

\[ \text{Beta}[a,b] \sim \theta^{a-1}\times(1-\theta)^{b-1} \]

So - if we want to create a uniform distribution for our prior, we can model this using a beta distribution by setting \(a\) and \(b\) to both equal to 1

\[ p(\theta) \sim \text{Beta}[1,1] \sim \theta^{1-1}\times(1-\theta)^{1-1}=\theta^{0}\times(1-\theta)^{0}= 1 \text{ (for all values of }\theta) \]

And we can model our likelihood using the binomial function

\[ p(X| \theta) = p(\text{H=}k| \theta) = \text{Binom}(n,p) \sim \theta^k\times(1-\theta)^{n-k} \]

And finally - the posterior is proportional to the prior multiplied by the likelihood:

\[ p(\theta) \times p(X| \theta) \sim [\theta^{a-1}\times(1-\theta)^{b-1}] \times [\theta^k\times(1-\theta)^{n-k}] = \theta^{a-1+k}\times(1-\theta)^{b-1+n-k} \]

Which we can purely parametize in a beta distribution (you can see that, the stronger the prior - and hence the larger a and b - the smaller the impact new evidence k and n has on the posterior distribution):

\[ \theta^{a-1+k}\times(1-\theta)^{b-1+n-k} \sim \text{Beta}[a+k,b+n-k] \]

(Note this is an approximation - the actual equality is specified later, as an example of conjugancy)

gamma_function <- function(n) {
  return(factorial(n-1))
}

beta_probability <- function(theta,a,b) {
  numerator <- ((theta)^(a-1))*((1-theta)^(b-1))
  denominator <- (gamma_function(a)*gamma_function(b))/gamma_function(a+b)
  return(numerator/denominator)
}

a <- 1
b <- 1

# uniform distribution
# prior <- dbeta(x = x, shape1 = 1, shape2 = 1)
x <- seq(0, 1, length=100+1)
prior <- sapply(x, FUN = function(theta) {beta_probability(theta,a,b)})

# binomial distribution
likelihood_f.binom <- sapply(x, FUN = function(theta) {binomial_probability(k,n,theta)})
likelihood_f.binom.int <- likelihood_f.binom
likelihood_f.binom.int[round(x*n,0) != x*n] <- NA

# beta distribution
likelihood_f.beta <- sapply(x, FUN = function(theta) {beta_probability(theta,k+1,n-k+1)*1/(n+1)})

posterior <- sapply(x, FUN = function(theta) {beta_probability(theta,a+k+1,b+n-k+1)*1/(n+1)})

plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = x, y = prior/sum(prior), name = 'Prior: Uniform = Beta(1,1)') %>%
  add_trace(x = x, y = likelihood_f.binom.int, name = 'Likelihood: ~ Binom(k=2,p=theta,n)', mode = 'markers') %>%
  add_trace(x = x, y = likelihood_f.beta, name = 'Likelihood: ~ Beta(p=theta,k,n-k)') %>%
  add_trace(x = x, y = posterior, name = 'Posterior: ~ Beta(a+k,b+n-k)') %>%
  layout(xaxis = list(title = 'Theta'))

As you can see, with a very ‘objective’ prior, the Bayesian posterior is very similar to the result derived from the frequentist evidence-only distribution (the Bayesian likelihood function).

Note here that this a good example of conjugancy - see the section below on what that means exactly.

In the two Bayesian examples above though, we either had a very weak prior (the uniform distribution) or a very strong prior (only two available values for theta, and 90% likelihood the coin is fair). Bayesian tends to most valuable where both the new data and past prior have more equal influence, and the posterior forming a distribution that sits somewhere between that estimated by the prior and likelihood function:

So - let’s say the last coin we tried flipping 10 times got heads 1.25 times (our prior) and the new coin flipped heads 2 times (our new evidence). We can then get the following equations:

\[ \text{Prior: } \text{Beta}[k_{1},n_{1}-k_{1}] = \text{Beta}[5,10-5] = \text{Beta}[5,5] \] \[ \text{Likelihood: } \text{Beta}[k_{2},n_{2}-k_{2}] = \text{Beta}[2,10-2] = \text{Beta}[2,8] \]

\[ \text{Posterior: } \text{Beta}[k_{1}+k_{2},n_{1}-k_{1}+n_{2}-k_{2}] = \text{Beta}[5 + 2,5 + 8] = \text{Beta}[7,13] \]

# prior <- sapply(x, FUN = function(theta) {beta_probability(theta,k1+1,n1-k1+1)/(n1+1)})
# likelihood_f <- sapply(x, FUN = function(theta) {beta_probability(theta,k2+1,n2-k2+1)/(n2+1)})
# posterior <- sapply(x, FUN = function(theta) {beta_probability(theta,k1+k2+2,n1-k1+1+n2-k2+1)/(n1+n2+2)})

# prior <- dbeta(x, k1+1, n1-k1+1)/(n1+1)
# likelihood_f <- dbeta(x, k2+1, n2-k2+1)/(n2+1)
# posterior <- dbeta(x, k1+k2+1, n1-k1+n2-k2+1)/(n1+n2+1)

prior <- dbeta(x, k1, n1-k1)
likelihood_f <- dbeta(x, k2, n2-k2)
posterior <- dbeta(x, k1+k2, n1-k1+n2-k2)

plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = x, y = prior, name = 'Prior belief: Horse') %>%
  add_trace(x = x, y = likelihood_f, name = 'Likelihood (New data): Donkey') %>%
  add_trace(x = x, y = posterior, name = 'Posterior Belief: Mule') %>%
  layout(title = "Updating Prior Beliefs")

The posterior curve is both taller and narrow than both the distributions estimated from the prior and the likelihood functions. This is because, by including both prior data and new evidence, our sample size has increased, so this reflects greater confidence that the observed number of heads,conditioned by the probability distribution of \(\theta\), will lie within a smaller interval.

Which leads on to Bayesian Credible intervals - the Bayesian equivalent to frequentist’s confidence intervals. Because Bayesian estimates distributions rather than fixed values, credible intervals are arguably more intuitive than confidence intervals: they reflect the a (95%) probability that of here the value generated by the random variable theta will fall within the range (rather than the frequentist interpretation that the true value will fall within the confidence interval 95% of the time).

Like a confidence interval, a credible interval is symmetric around the mean, and spans the x-axis enough to cover 95% of the area under the probability distribution curve. And that is good if the posterior distribution is also symmetric - but as can be seen in our example above, that is not necessarily the case, where the posterior is skewed.

If the distribution is skewed then, it could be better to compute the region of highest posterior density (HPD) - which by definition will be the narrowest interval across the potential values for theta, but cover 95% of the area under the probability curve. In other words, we still maintain the same 95% probability of theta lying within the region, but the region is the narrowest it possibly can be.

So to do this, we want to find the interval between \(\theta_{low}\) and \(\theta_{high}\) such at:

\[ \int_{\theta_{low}}^{\theta_{high}} p(\theta|x) d\theta= 1-\alpha \]

But finding this can be tricky, particularly with computing closed intervals. The code below shows a brute force solution. The way to think about this is as follows: * Credible intervals: Given the mean of the distribution, expand the range outwards equally until 95% of the area under the curve is captured. * Highest posterior density: Imagine a horizontal line moving slowly towards the x axis, intercepting the distribution twice. Once the

(NOTE THE GRAPH BELOW IS WRONG, NEEDS FIXING)

find_binomial_credible_interval <- function(k, n, p, learning_rate = 0.5, initial_width = 0.1, criterion = 0.1, max_iterations = 1000, precision = 4) {
  bit <- initial_width/2 # amount to push away from mean
  for(i in 1:max_iterations) {
    amount <- pbeta(p+bit, k+1, n-k+1) - pbeta(p-bit, k+1, n-k+1) # sum the auc between p-bit and p+bit
    bit <- bit - learning_rate*(amount - (1-criterion)) # decrease if auc > 0.95 or increase if auc < 0.85 
    if(abs(amount - (1-criterion)) < 1/(10^precision)) {
      break()
    }
  }
  return(c(p-bit,p+bit))
}

find_binomial_highest_start_point <- function(precision = 4,k,n) {
  x <- seq(0,1,10^-precision)
  pdf <- dbeta(x,k,n-k)
  return(x[max(pdf) == pdf])
}

find_binomial_hpd <- function(k,n,criterion=0.1, initial_movement_size = 0.1, max_iterations=10^3,learning_rate =  0.005, precision = 4) {
  p <- find_binomial_highest_start_point(precision,k,n)
  left_loc <- p
  right_loc <- p
  bit <- initial_movement_size
  if(residual > 0) { # need to increase width
    if(pbeta(left_loc + bit, k+1, n-k+1) > pbeta(right_loc + bit, k+1, n-k+1)) { # left side pdf is higher
      left_loc <- left_loc-bit # move left
    } else { # right side pdf is higher
      right_loc <- right_loc+bit # move right
    }
  } else {
    if(pbeta(left_loc + bit, k+1, n-k+1) > pbeta(right_loc + bit, k+1, n-k+1)) { # left side pdf is higher
      right_loc <- right_loc+bit # move right
    } else { # right side pdf is higher
      left_loc <- left_loc-bit # move left
    }
  }
  residual <- (1-criterion) - (pbeta(right_loc, k+1, n-k+1)-pbeta(left_loc, k+1, n-k+1)) # difference to e.g. 0.95
    if(abs(residual) < 10^-precision) {
      break()
    } else {
      bit <- bit - learning_rate*residual
    }
  return(c(left_loc, right_loc))
}

credible_intervals <- find_binomial_credible_interval(k, n, p, max_iterations = 10^3, criterion = 0.1)
highest_posterior_density <- find_binomial_hpd(k,n, max_iterations = 10^3, criterion = 0.1)

Error in find_binomial_hpd(k, n, max_iterations = 10^3, criterion = 0.1) : 
  object 'residual' not found

Some Bayesians argue though that credible intervals and HPDs are only really something that are useful to compare to frequentist interpretations, as the distribution is already given by the posterior.

Conjugancy (beta-binomial)

Now to get geeky - you might have already noticed how eerily similar the beta distribution is to the binomial distribution from the formula - and in fact you can use the beta distribution as a continous function equivalent to the binomial function.

\[ p(\text{H=}k| \theta) \sim \theta^k\times(1-\theta)^{n-k} \sim \text{Beta}[k+1,n-k+1] \]

Both binomial and beta likelihoods are scaled though (normalized), so that the sum of their probabilites total 1. For the beta distribution, the scaling is achieved like this:

\[ Beta[\alpha,\beta] = \theta^{\alpha-1}\times(1-\theta)^{\beta-1} \times \frac{1}{\text{B}(\alpha,\beta)} \] Where

\[ \frac{1}{\text{B}(\alpha,\beta)} = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \times \Gamma(\beta)} = \frac{(\alpha+\beta-1)!}{(\alpha-1)! \times (\beta-1)!} \]

So if successes \(\alpha = k+1\) and failures \(\beta = n - k + 1\), then:

\[ \frac{(\alpha+\beta-1)!}{(\alpha-1)! \times (\beta-1)!}=\frac{((k+1)+(n-k+1)-1)!}{((k+1)-1)! \times ((n-k+1)-1)!}=\frac{(n+1)!}{(k)!\times(n-k)!} = {n+1 \choose k } = (n+1) \times {n \choose k } \]

Now we define the binomial distribution for comparison:

\[ \text{Binom}(n,p) = {n \choose k } [\theta^{k}\times(1-\theta)^{n-k}] \]

And thus:

\[ Beta[k+1,n-k+1] = (n+1) \times \text{Binom}(n,p) \]

This relationship helps us combine binomial likelihoods and beta priors to create a posterior in what is called conjugancy. Conjugancy occurs when the posterior distribution is in the same family of probability density functions as the prior belief, but with new parameter values (updated to reflect what has been learned from the data). The posterior comes from the same family as the prior, rather than the data’s distribution.

All the three coin examples above were examples of conjugancy. For example, the second had a prior that was uniform, expressed as a beta distribution, whereas the likelihood function was a binomial distribution. The posterior distribution was thus beta, since it is the same family of probability density functions as the prior belief (beta), but with new parameter values learnt from the probability mass function derived from the new data. Note that the posterior did not follow a binomial distribution, but followed a beta one, since this is what the prior used. Conjugates are thus very useful in helping derive the posterior distribution.

Without conjugancy, one has to do the integral, which is often impossible to actually evaluate, and is one of the reasons bayesian was not popular in the 20th century (before computing allowed researchers to compute integrals numerically). Nowadays, sampling methods such as Markov-chain monte-carlo (MCMC) are often used to simulate the posterior distribution.

Conjugancy (gamma-poisson)

As well as the beta-binomial conjugate family, we also have the gamma-poisson conjugate family.

The poisson function is specified as:

\[ P(X=k)=\frac{\lambda^{k}}{k!}e^{-\lambda} \text{ for } k=0,1,..., \]

Whereas the gamma family is related to the gamma function (witk \(k,\theta\) and \(\alpha,\beta\) parameterizations listed below):

\[ f(x)=\frac{1}{\Gamma(k)\theta^k}x^{k-1}e^{-x/\theta}=\frac{\beta^{\alpha}}{\Gamma(k)}x^{\alpha-1}e^{-\beta x} \] And hence we can see that that where \(n\) is a positive integer, and \(\theta=1\), we observe the special case that \(\Gamma(k)=(k-1)!\) so that

NEED TO WORK ON THIS BIT

\[ \frac{1}{\Gamma(k)\times\theta^k}x^{k-1}e^{-x/\theta}= \frac{1}{(k-1)!\times1^k}x^{k-1}e^{-x/\theta}= \frac{\lambda^{k}}{k!}e^{-\lambda} \]

Thinking about Bayesian in terms of predicting things

There are two sources of uncertainty when building statistical models to predict things:

• Uncertainty due to the fact any future value is itself a random event, \(P(y|\theta)\) (because future events are essentially samples that have error conditioned by \(\theta\))
• Uncertainty in the parameter values which have been estimated on the basis of past data, \(P(\theta|X)\) (for example, a coefficient estimated in a regression model)

Classical models deal with point 1, making a point estimate in the future based on the ‘best’ parameters estimated from past data. This is equivalent to forming a likelihood (i.e. incorporating new data). What they don’t do though is think about point 2, that the seemingly optimum model estimated itself might be incorrect since the true values of \(\theta\) are not fixed, and this is equivalent to not building in any prior beliefs. In a way, by giving a point estimate, the models generate a false sense of precision, and the confidence intervals arounds any estimate is only predicated on the idea there is sampling randomness, rather than randomness associated with \(\theta\) itself.

In other words, classical models seek to best fit the new evidence (traditionally through maximum likelihood estimation for generalized linear models), whereas bayesian models incorporate prior beliefs as well.

Imagine then, if we want to estimate the impact of x on y to make a future prediction, we can do the following:

\[ f(y|x)= \int f(y|\theta,x)f(\theta|x) \partial(\theta) \]

Where \(f(y|\theta,x)\) represents the modelling random sampling error as in frequentist regression (the likelihood function), and \(f(\theta|x)\) representing the error arising from \(\theta\) being a random variable rather than being treated as a fixed quantity.

Now why is it an integral? Let’s think of the coin example above where there are only two types of coins in the economy. If we want to know what is the probability we get heads, rather than the probability that we have a fair or a biased coin.

\[ P(\text{heads}|k=2) \\ = P(\text{heads}|\text{fair},k=2) \times P(\text{fair}|k=2) + P(\text{heads}|\text{biased},k=2) \times P(\text{biased}|k=2) \]

As before, we summarise the probability mass function relating to \(\theta\) as being 0.5 (fair) for 90% of coins in the economy, and $0.2 (biased) for 10% of these coins, which gives us the following:

\[ P(\text{heads|k=2}) \\ = \sum_{\theta}{P(\text{heads}|\theta\text{,k=2}) \times P(\theta|\text{k=2}) } \]

And by taking into account all the possible values that theta can take, we are essentially integrating along the

\[ P(\text{heads|k=2}) = \int P(\text{heads}|\theta,k=2) \text{ }P(\theta|k=2) \text{ } \partial(\theta) \\ \text{where }\theta = \begin{cases} 0.2 & 10\% \text{ of the time}\\ 0.5 & 90\% \text{ of the time}\\ \end{cases} \]

NEED TO RE-WRITE EVERYTHING BELOW THIS

!!!!!!!!!!!!! !!!!!!!!!!!!! !!!!!!!!!!!!! !!!!!!!!!!!!! !!!!!!!!!!!!! !!!!!!!!!!!!! !!!!!!!!!!!!! !!!!!!!!!!!!! !!!!!!!!!!!!!

Now because this equation is done through conjugancy, this is easy to evaluate. Another way then of looking at it for this example is by rearranging bayes theroem:

\[ P(y|X) = \frac{P(y|\theta,X) \times P(\theta|X)}{P(\theta|y,X)} = \frac{P(\text{heads}|\text{fair},k=2) \times P(\theta|X)}{P(\theta|y,X)} \]

$$ P(|X) = = \

$$

And hence it can be seen that to get the full probability distribution associated with getting heads, you need to assess all possible values that theta could take.

\[ P(y|X) = \frac{P(y|\theta,X) \times P(\theta|X)}{P(\theta|y,X)} = \frac{P(\text{heads}|\theta,k=2) \times P(\theta|k=2)}{P(\theta|\text{heads,k=2})} \]

$$ P() = + \

$$

As before, we summarise the probability mass function relating to \(\theta\) as being 0.5 (fair) for 90% of coins in the economy, and 2 (biased) for 10% of these coins, which gives us the following:

\[ P(\text{heads|k=2}) = P(\text{heads}|\theta\text{,k=2}) \times P(\theta|\text{k=2}) \]

\[ P(\text{heads}) = P(\text{heads}|\text{fair}) \times P(\text{fair}) + P(\text{heads}|\text{biased}) \times P(\text{biased}) \]

\[ P(X,\theta|k=2) = \frac{P(\text{heads}|\theta,k=2) \times P(\theta|X)}{P(heads|k=2)} \]

\[ P(y|X) = \frac{P(y|\theta,X) \times P(\theta|X)}{P(\theta|y,X)} \]

Let’s say we believe the sampling error of y to follow a normal distribution, conditioned by \(\theta\) (capturing the first source of uncertainty), which we would normally solve through MLE:

\[ f(y|\theta) = \frac{1}{\sigma\sqrt{2\pi}}\text{exp}(\frac{(y-\theta)^2}{2\sigma^2}) \] \[ f(y|\theta,X) = \frac{1}{\sigma\sqrt{2\pi}}\text{exp}(\frac{(y-X\beta)^2}{2\sigma^2}) \]

Then we just need to derive \(f(\theta|x)\) in order to make a future prediction through bayes (and this term captures the error associated with our estimation of theta from the exogenous variables).

Now for frequentists, the pdf \(f(\theta|x)\) is equal to 1 at the true value of theta, and zero elsewhere, as theta is fixed (so there is only one perfect value).

Though for Bayesians, the pdf is conditioned by prior beliefs about \(\theta\). These are represented as a posterior approximated from the product of the likelihood of observing the exogenous variables given \(\theta\), and the prior beliefs associated with \(\theta\):

\[ f(\theta|x) \sim f(\theta) \times f(x|\theta) \] To capture the randomness of the \(\theta\), we might also believe our prior \(f(\theta)\) follows a normal distribution, but with mean \(E(X)=b\) and variance \(Var(X)=d^2\). As such, the distribution thus described is \(X \sim N(b,d^2)\)

Then the prior can be derived from the normal pdf \(x\):

\[ f(\theta) = \frac{1}{d\sqrt{2\pi}}\text{exp}(\frac{(\theta-b)^2}{2d^2}) \]

Now let’s take the observed data \(x=(x_1,x_3,...,x_n)\) and treat it as a random sample of size \(n\) of a random variable \(X\) with mean \(E(X)=\theta\) and variance \(Var(X)=\sigma^2\), and follows a normal distribution thus described as \(X \sim N(\theta,\sigma^2)\)

Then the likelihood for every obervation \(x_i\) from the whole sample \(x\):

\[ f(x|\theta) = \prod_{i=1}^{n} f(x_i|\theta) = \prod_{i=1}^{n} \frac{1}{\sigma\sqrt{2\pi}}exp(\frac{(x_i-\theta)^2}{2\sigma^2}) \]